'''
Created on Dec 9, 2010

@author: oabalbin
'''
import exome.gatk_cluster.cluster_jobs_header as jh


def snps_calls_VarScan(list_bam_normals,list_bam_tumors, snps_calls_file, use_mem, path_to_varscan):
    '''
    java -jar VarScan.v2.2.jar somatic normal.pileup tumor.pileup output.basename
    The above command will report germline, somatic, and LOH events at positions 
    where both normal and tumor samples have sufficient coverage (default: 8). 
    If the --validation flag is set to 1, it also report non-variant (reference) 
    calls at all positions with sufficient coverage, for the purpose of refuting 
    false positives with validation data. 
    
    Writes:
    output.snp
    output.indel
    '''
    
    varscan_command=path_to_varscan+'VarScan.v2.2.3.jar'
    args = ['java', '-Xmx'+str(use_mem)+'m', '-jar', varscan_command, 'somatic', 
            '--validation', '1']+ list_bam_normals + list_bam_tumors + [snps_calls_file]
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command


def process_somatic_calls(snps_calls_file, path_to_varscan):
    '''
    java -jar VarScan.v2.2.jar processSomatic [output.snp]

    Writes:
    output.snp.Somatic.hc     (high-confidence Somatic mutations)
    output.snp.Somatic.lc     (low-confidence Somatic mutations)
    output.snp.Germline     (sites called Germline)
    output.snp.LOH         (sites called loss-of-heterozygosity, or LOH)    -
    '''
    varscan_command=path_to_varscan+'VarScan.v2.2.3.jar'
    args = ['java', '-jar', varscan_command, 'processSomatic', snps_calls_file]
    
    args= [a.replace(',',';') for a in args]
    command = ",".join(args).replace(',',' ').replace(';',',')

    return command

    
    
def main_VarScan_caller(cohort_dict, run_dict):
    '''
    '''
    
    path_to_varscan=run_dict['path_to_varscan']
    short_wtime=run_dict['short_wtime']
    extra_mem=run_dict['use_mem']
    out_dir = run_dict['out_dir']
    node_processors = run_dict['num_cores']
    single_processor=1
    my_email = run_dict['myemail']
    
    jobn='svs'
    i=0
    for cohort, cohort_samples in cohort_dict.iteritems():
        # Note that in this case there list_bam_files_normal
        # correspond to a merged bam file of normals or tumors
        # for that cohort
        jobn=jobn+str(i)
        list_bam_files_normal = cohort_samples['normals']
        list_bam_files_tumor = cohort_samples['tumors']
        
        
        snps_calls_file = out_dir+cohort+'.varscan_snps'
        command = snps_calls_VarScan(list_bam_files_normal,list_bam_files_tumor, snps_calls_file, extra_mem,
                                      path_to_varscan)
        
        jobidvc = jh.qsub(jobn+'_vc', command, single_processor, cwd=None, walltime=short_wtime, pmem=extra_mem, 
                           deps=None, stdout=None, email_addresses=my_email)
        
        command = process_somatic_calls(snps_calls_file+'.snp', path_to_varscan)
        
        jobidpvc = jh.qsub(jobn+'_pvc', command, single_processor, cwd=None, walltime=short_wtime, pmem=extra_mem, 
                           deps=jobidvc, stdout=None, email_addresses=my_email)
        
        # Look how to convert the final varscan output into a vcf format
        
        i+=1
        
        return jobidpvc
        
if __name__ == '__main__':
             
    cohort_dict = {'Prost':{'normals':['/nobackup/med-mctp/oabalbin/test/Prost.normals.sam_snps.mpileup.vcf'],
                            'tumors':['/nobackup/med-mctp/oabalbin/test/Prost.tumors.sam_snps.mpileup.vcf']}}
    
    run_dict = {'path_to_gatk':'/nobackup/med-mctp/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/',
                 'resources_folder':'/nobackup/med-mctp/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/resources/', 
                 'rscipt_path':'/home/software/rhel5/R/2.10.1-gcc/bin/Rscript',
                 'path_to_picard':'/nobackup/med-mctp/sw/bioinfo/picard/picard-tools-1.35/', 
                 'path_to_bwa':'/nobackup/med-mctp/sw/bioinfo/alignment/bwa/bwa-0.5.8a/',
                 'path_to_sam':'/nobackup/med-mctp/sw/bioinfo/samtools/samtools-0.1.10/',
                 'path_to_vcftools':'/nobackup/med-mctp/sw/bioinfo/vcftools/',
                 'path_to_varscan':'/nobackup/med-mctp/sw/bioinfo/varscan/',
                 'use_mem':8000, 'num_cores':6,
                 'nmismatches':2,
                 'short_wtime':'24:00:00',
                 'long_wtime':'60:00:00',
                 'ref_genome':'/nobackup/med-mctp/sw/alignment_indexes/gatk/hg19/hg19.fa',
                 'ref_genome_bwa':'/nobackup/med-mctp/sw/alignment_indexes/bwa/hg19/hg19.fa',
                 'snpdb_file':'/nobackup/med-mctp/sw/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed.vcf',
                 'indeldb_file':'/nobackup/med-mctp/sw/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed.vcf',
                 'target_exons':'/nobackup/med-mctp/sw/alignment_indexes/agilent/hg19/exome/Agilent_SureSelect_All_Exon_G3362.v2.hg19.bed',
                 'recal_analysis_outputdir':'/nobackup/med-mctp/oabalbin/test/recal_analysis/',
                 'temp_dir':'/nobackup/med-mctp/oabalbin/test/temp/',
                 'qsubfile':'/nobackup/med-mctp/oabalbin/test/',
                 'out_dir':'/nobackup/med-mctp/oabalbin/test/',
                 'myemail':['alebalbin@gmail.com']
                 }
    
    main_VarScan_caller(cohort_dict, run_dict)



        

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
        
